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Abstract. Computation of high frequency solutions to wave equations is impor- 
tant in many apphcations, and notoriously difficult in resolving wave oscillations. 
Gaussian beams are asymptotically valid high frequency solutions concentrated 
on a single curve through the physical domain, and superposition of Gaussian 
beams provides a powerful tool to generate more general high frequency solutions 
to PDEs. An alternative way to compute Gaussian beam components such as 
phase, amplitude and Hessian of the phase, is to capture them in phase space by 
solving Liouville type equations on uniform grids. Following j^j we present a sys- 
tematic construction of asymptotic high frequency wave fields from computations 
in phase space for acoustic wave equations; the superposition of phase space based 
Gaussian beams over two moving domains is shown necessary. Moreover, we prove 
that the fc-th order Gaussian beam superposition converges to the original wave 
field in the energy norm, at the rate of t^+^r- in dimension n. 



Contents 

1. Introduction 1 

2. Phase space based Gaussian beam Ansatz 3 

3. Recovery of the high frequency wave fields 5 

4. Control of initial error l7 

5. Propagation of the approximation error 10 

6. An example 13 

7. Higher order Approximations 14 
Acknowledgments 17 
References 17 



1. Introduction 

This is the continuation of our project, initiated in [3], of developing a rigorous 
recovery theory for high frequency wave fields from phase space based computations. 
Here we focus on the wave equation 

(1.1) Pu := [d"^ - c{xyA]u = 0, (x, t) E M" x R, 

where c{x) is a positive smooth function, with highly oscillatory initial data 

(1.2) u{x, 0) = Ai^{x, e)e^^'"(^)/S Ut{x, 0) = Bi^{x, e)e^^-(^)/^ 
Date: June 15, 2009. 
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The initial phase 5;^ G C°°(M"), and the amphtudes v4in,An e C(5^(M") have the 
following asymptotic expansions: 

(1.3) An : = At\x) + eA«(x) + e^J^f (x) + ■ • • , 

(1.4) i?in : = e-'Bl;'\x) + i?J°)(x) + ei?«(x) + ■ • • . 

The small parameter e represents the typical wave length of oscillations of the ini- 
tial data. Propagation of oscillations of wave length e causes mathematical and 
numerical challenges in solving high frequency wave propagation problems. 

In this article we are interested in the construction of globally valid asymptotic 
wave fields and the analysis of their convergence to the true solutions of the initial 
value problem. A general discussion of this problem and background references are 
given in the introduction to |3j. We have two objectives: 

i) to present the construction of asymptotic solutions as superpositions over 
phase space; 

ii) to estimate the difference between the exact wave fields and the asymptotic 
ones. 

The construction for (i) is based on Gaussian beams (GB) in physical space con- 
structed similarly to those given for wave equations in [B] , but here the construction 
is carried out by solving inhomogeneous Liouville equations in phase space. While 
the result is no longer a superposition of asymptotic solutions to the wave equa- 
tion (1.1), the superposition is nonetheless asymptotic. We consider superpositions 
over two suhdomains moving with two Hamiltonian flows, respectively, and show 
that they are asymptotic solutions by relating them to the Lagrangian superposition 
through two time-dependent symplectic changes of variables. An argument of this 
type was used for the Helmholtz equation in [2J. 

For (ii), as in [6j, we use the well-posedness theory for (1.1), i.e. the continuous 
dependence of solutions of Pip = / on their initial data and /. Thus, the sources of 
error in the Gaussian beam superposition for the initial value problem are the error 
in approximating the initial data and the error in solving the PDE. There are some 
differences between the acoustic wave equation and the Schrodinger wave equation. 
For example, the caustics that can form are weaker. 

In summary, our phase space based Gaussian beam superposition is expressed as 

/ '^PGB 

Jn+{t) Jn~{t) 

where X = {x,p) denotes variables in phase space M^*^, Q{0) is the domain where we 
construct initial Gaussian beams from the given data, and f2^(t) is the image of fl{0) 
under the Hamiltonian flow for H{x,p) = ±c{x)\p\. The functions u^Qg{t^ y, X) are 
constructed using the phase space based Gaussian beam Ansatz, and Z{n, e) ~ g""/^ 
is a normalization parameter. Our result shows that for the k—th order phase space 
Gaussian beam superposition, the following estimate holds 

(1.6) Uu^ - u){t, ■)|U < 11(^^(0, ■) - uU-)\\e + |^^(0)|e^+V, 

where ||e|||; := y /]Rn[c^^|etp + iVi^e^jda;. Here and in what follows we use A < B 
to denote the estimate A < C5 for a constant C which is independent of e. 



(1.5) u'{t,y) = Z{n,e) 
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;i.7) u\t,y) = Z{n,t) 



For the initial data of the form (Ain(x, e), -Bin(a:, e))e*'^'"*^'^''/^ we need a superposi- 
tion over an n-dimensional submanifold of phase space. The asymptotic solution is 
then represented as 

/ UpQQ5{w~^)dX + / Up(jg6{w^)dX , 
Jn+{t) JQ^it) J 

where is obtained from the Liouville equation 

dtw + Hp ■ V^w - ■ Vpw = 0, w{0, X) = p - V^Si^ix), 
with H{x,p) = ±c(x)|p|, respectively. Our result shows that 

(1-8) ll(«^-«)(t,-)IU<et+^. 

Here the exponent k/2 reflects the accuracy of the Gaussian beam in solving the 
PDE. It will increase when one uses more accurate beams. The exponent indi- 
cates the damage done by the caustics. 

We now conclude this section by outlining the rest of this paper: in Section 2 we 
start with Gaussian beam solutions in physical space, and define the phase space 
based GB Ansatz through the Hamiltonian map. Section 3 is devoted to a recovery 
scheme through superpositions over two moving domains. The total error is shown 
to be bounded by an initial error and the evolution error of order 6(3-n)/4 Qgntrol of 
initial error is discussed in Section 4. Convergence rates are obtained for first order 
GB solutions in Section 5. In Section 6 we present an example to illustrate these 
constructions. Extensions to higher order GB approximations are given in Section 
7. 

2. Phase space based Gaussian beam Ansatz 

As is well known, the idea underlying Gaussian beams [5J is to build asymptotic 
solutions concentrated on a single ray path in Mt x M". This means that, given a 
ray path 7 parameterized by {t, x{t)), one makes the ansatz 

(2.1) u^{t,y)=A{t,y,e)e'''^''yy% 

where $(t, x{t)) is real, and /m{$(t, y)} > for ?/ 7^ x(t). The amplitude is allowed 
to be complex and has an asymptotic expansion in terms of e: 

A{t, y, e) = Ao{t, y) + eA^it, y) + ■ ■ ■ + e'^v4^(t, y). 

We wish to build asymptotic solutions to Pu(t,y) = 0, i.e., we want Pu'' = O(e^). 
Substituting from (12. ip and grouping terms multiplied by the same power of e, we 
obtain the equations of geometric optics: 

(2.2) P[A{t, y, e)e^*(*'^)A] = (^f^ c,{t, y)e^^ e^*(*'^)/% 

where for G{t,y) = - c^|Vy$p, 

c_2(t,y) = -G{t,y)Ao, 
c.i{t,y) = 2iLAo + G{t,y)Ai, 

ci_,{t, y) = 2iLAi + G(t, y) A+i + P[Ai-il / = 1, ■ ■ ■ A - 1. 
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Here L is the linear differential operator, 

Since e**/^ decays rapidly away from 7, to make P(Ae**/^) = 0(e^'^) for a given 
M G Z, we only need to make Cj vanish on 7 to sufficiently high order. In this 
work we discuss mainly the lowest order Gaussian beam solutions, followed by an 
extension to higher order Gaussian beam superpositions in Section 7. 
We begin with c_2 = 0, i.e. (7 = 0. This leads to two eikonal equations 

(2.3) (9t<l> + V^*) = 0, H{x,p) = ±c{x)\p\. 

The leading amplitude solves 

(2^4) a,^ + ff,.V.A = ^^^lM_^ 

We continue to denote the phase space variable as X = and let Xq = (xo,Po) 

denote the initial state. Then the equations for the bicharacteristics X = X^{t, Xq) 
originating from Xq at t = are 

(2.5) ^^^^'^0) = V{X{t,Xo)), X(0,Xo) = Xo. 

The vector field V = {Hp, —Hx) is divergence free, and hence this flow preserves the 
volume on phase space. 

From now on we include the initial data Xq as a parameter in the phase: $ = 
$(t,l/;Xo) and the amplitude: A = A{t,y;XQ). We apply Taylor expansion of the 
phase $ and the amplitude A about x = x(t,Xo) to obtain 
(2.6) 

$(t, y; Xo) = Sit; Xo)+p(t, Xq) (y-x(t, Xo)) + ^(y-x(t, Xo))^M(t; Xo)(y-x(t, Xq)), 

with p{t, Xq) = dy^{t, x{t, Xo); Xq) and 

S{t; Xo) = <l>(t, x{t, Xo); Xo), M(t; Xo) = ^^^(t, x{t, Xo); Xo). 

For the amplitude we set A(t,y;Xo) = A{t;Xo) with A(t;Xo) = A(t,x(t,Xo);Xo). 
Then we get the equations along the curve 7 for S 

(2.7) ^5(t;Xo) = 0, S{0;Xo) = S^xo), 

at 

and the Hessian M 

(2.8) j^M{t; Xo) + + H^pM + MHp, + MHppM = 0, M(0; Xo) = Min(xo). 

Using the eikonal equation (9^$ + H{x, V$) = twice, we see that 

P[$] = dt[-H{x, V$)] - c^A* = Hp-Hx + HpMHp - c^Tr{M). 
This with ( 12. 4p shows that the amplitude along the ray, A(t;Xo), satisfies 

(2.9) j^A{t; Xo) = ^[Hp-H, + HpMHp - c^Tr{M)] , A{0; Xo) = A,(a;o). 

We have introduced this form of the transport equation because it is easier to trans- 
late to Eulerian coordinates. The essential idea behind the Gaussian beam method 
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is to choose some complex Hessian Mjn initially so that M remains bounded for 
all time, and its imaginary part is positive definite. Equation (12.91) shows that the 
amplitude A{t] Xq) will also remain bounded for all time. 

The above construction ensures that the following GB Ansatz is an approximate 
solution 

UGsit, y; Xo) = u%^{t, y; Xq) + u:^j^{t, y; Xq), 

where 

u^^{t, y- Xo) = A^it- Xo) exp {^-^^H. V\ Xq] 

where both A^{t] Xq) and y; Xq) are computed from (12.91) and (12.61) with H = 
±c{x)\p\, respectively. 

Here ^^(0; X) are to be chosen so that a superposition will match the initial data 

(M,Mi)|t=o= (An,Si,)e^^-/^ 

to leading order. For this matching we need (for X = (x, V5'in(x))) 

A^iO;X) + A-iO;X) = A^ix), 

'-A+iO; X)dt^+{0, X- X) + '-A' (0; X)dt^- (0, x; X) = -^b[-^^ (x). 

In the second relation we took only the leading term in e~*'^'"/^Mt. Since the two 
Hamiltonians have different signs, 

$±(0,x;X) = 5in(x) and 9t$±(0,x;X) = Tc(x)|V,.^in(x)|, 

the second relation gives 

(2.10) A+(0;X) - A-(0;X) = , '^j"^,^ .. . 

c(x)|Va;Sin(x)| 

Hence solving for A^ we have 

1 / .(0). . . iBt^Hx) 



[2.11) ^^(0;X) = -M[„^^(x)± 



cix)W :rS\Jx) 



Note that we could simplify the superposition by taking some special initial data 
such that B^~^\x) = — zA[°-'(x)c(x)|V5'in(x)|. The advantage of these special choices 
is that we do not need a sum of two Gaussians to approximate the solution. We also 
note that for given initial of order 0(1), i.e., B^ = 0, we see that A^{0;X) = 



1 /l(O) 



2 An (^)- 



3. Recovery of the high frequency wave fields 

Since the wave equation we consider is linear, the high frequency wave field u at 
(t, y) in physical space is expected to be generated by a superposition of neighboring 
Gaussian beams 

(3.1) u'{t,y) = Zin,e) [ UGBit,y; Xo)dXo, 

Jn(o) 
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where Q{0) is a bounded open set containing 

{Xo : xo G supp(Ain) U supp(5in), po e Tange{d^Sin)}. 

The normahzation parameter Z(n, e) ~ e""/^ is determined by matching initial data 
against the Gaussian profile. 

Since the flows Xq) are volume preserving in phase space, 

Using X = X^(t,Xo) and their inverses Xq = XQ(t,X), we obtain our Gaussian 
beam Ansatz in phase space 

From (13. ip it follows that 



Z{n,e) / [u^j^{t,y;Xo) +UQB{t,y]Xo)]dXo 
Jn(o) 



(3.2) = Z{n, 

where 



/ UpQ^{t,y,X)dX + / Upfj^ {t,y,X)dX 
Jn+{t) Jn-{t) 



n^(t) = x^{t,n{o)). 

Each phase space Gaussian beam has the form 

(3.3) upGBit, y, X) = Ait, X) exp (^^^t, y, X)^ , 
where 

(3.4) $(t, y, X) = S{t, X)+p-{y~x) + ^{y- xf M{t, X){y - x). 

Note that though UpQp{t, y, X) are no longer asymptotic solutions of the wave equa- 
tion in {t,y), their superpositions over the moving domains f2^()f:) in X remain as- 
ymptotic solutions. 

Let C be the Liouville operator defined by 

(3.5) C:=dt + V-Vx. 

If w{t, X) is the phase space representative of w(t; Xq) in the sense that w{t; Xq) = 
w{t, X(t, Xq)) for any t > 0, then 

jwit;XQ) = Cw{t,X). 

Hence from the Lagrangian formulation of equations for {S, M, A) we obtain PDEs 
for (5,M,i) in (ETD, (EED and (EH): 

(3.6) C{S) = 0, S{0,X) = SUx), 

(3.7) C{M) + H^^ + H^pM + MHp^ + MHppM = 0, M(0, X) = Min(x), 

(3.8) C{A) = -^\Hp-H^ + HpMHp - c^Tr{M)] , i(0,X) = A,(x), 
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where H{x,p) = c{x)\p\ or H{x,p) = —c{x)\p\. The heart of the matter is equation 
(13. 7p . It is known from [1] that, if M^^ is symmetric and the imaginary part of M^^ 
is positive definite, then a global solution M to (13. 7p is guaranteed and has the 
properties: 

i) M = M^, and 

ii) Im{M) is positive definite for all t > 0. 

There are several ways of computing M. Following [1] (see also [31 Section 7]), we 
use a level set method to construct the Hessian: 

(3.9) M = -g,{gpr\ 

where g = (f)i{t,X) + i(j)2{t,X) with (pi obtained by solving the Liouville equation 

C{<P) = 0. 

From the well-posedness theory of the wave equation we have the following. 

Lemma 3.1. Let u satisfy P[u\ = in [0, T] x M" with {u,Ut) given att = 0, and 
let be an asymptotic solution. Then the error e = u'^ — u satisfies 

(3.10) ||e(t)|U< ||e(0)|U + e / \\c-'P[u^]\\^,dT, 

Jo 

where \\e\\E = and 

E.= '^[ [c-'\e,\'+\V.e\']dx. 

Proof. Since we start with the data with compact support, at any finite time the 
support of the solution remains bounded (due to finite speed of propagation for the 
wave equation). 

Let e = u"^ — u. Then from P[u\ = 

P[e] = P[u']-P[u] = P[u']. 

We now have 

d f 

-E{t) = eM [c^^etctt + Ve ■ Ve^] dx 

at Jjln 

= e^ f [V ■ (e^Ve) + c-^etP[u']] dx 

< \\c-^et\\^, \\c-^P[u%\^, < eV2E\\c-^P[u%^, . 
This, upon integration in time, leads to the desired estimate. □ 

4. Control of initial error 

For the initial phase S'in, we set po = ^xSin^xo) and form the Lagrangian super- 
positions 

u'{t,y) = Z{n,e) / UGB{t,y; Xo)S{po - V^Si,,{xo))dXo. 
Jn{o) 
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In order to track the deformation of the surface p — VxSin{x) = as time evolves, 
we introduce two level set functions w = w^(t, X) such that 

(4.1) C[w] = 0, wiO,X)=p-V,SUx), 

with H{x,p) = ±c{x)\p\. Here gives 02 needed in (13. 9 p and 0i can be obtained 
from solving the respective Liouville equation with 0i(O,X) = x. 

Using the volume preserving maps X = X'^{t,XQ), leads to the Gaussian beam 
superposition in phase space 



/ UpQQ5{w~^)dX + / Upqq6{w )dX 
Jn+a) Jn-(t) 



(4.2) u'{t,y) = Z{n,e) 

-(*) Jn-it) 

where Q'^lt) = X^{t,Q{0)). Our choice of initial data for the beams in this super- 
position will be made to match the initial data in (11.21) . Set 

(4.3) /(0) = {x: {x,p)eQ{0), p = V,SM}. 

We now use the Lagrangian formulation of the GB superposition to match the initial 
data. 



(4.4) u^{t,y) = Z{n,e) UGBit,y; Xo)dxo. 

Jm 

Here and in what follows we use ucBit, y, xq) for UGsit, y, xo, VxSin{xo)). If we take 
S^{0;xo) = Si^ixo), M^{0;xo) = dlS^J^Xo) + if3I with /5 > as well as A±(0;Xo) 
as defined in (12. lip , then 

u^{0,y) = Z{n,e) [ A;°)(xo)e^*(°'^^^«)/^rfxo, 
Ji(o) 



where 



<!>{0,y-xo) = T^'[S,^]{y) - ^\y - xo\'. 



Here Tj^[S]{y) denotes the j^^ order Taylor polynomial of 5* about x at the point y. 
Setting 

Z{n,e) 

we have 



27re 



u'{0, y) = / A?(^o)e^[^"[^-l(^)]ir xo-y,4^] dxo, 



7(0) 



2P 



.|2 



where K{x, r) = ^^^^^^^^ e 4^ is the usual heat kernel, satisfying limitT-|o-f^(a;, t) 
6{x) as distributions on M", and 



K{x -y,T)dx = 1, Vr>0, yGR". 

J X 



On the other hand the initial wave field is 



u{0, y) = AL(y)e^^-(^)/^ = AL(y)e^^-(^)/^K (^x - y, dx. 
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Both the phase and amphtude in the integrand can be approximated by their Taylor 
expansion when |x — ?/| is small, say |x — ?/| < e^^^, and the integral over the com- 
plement of this neighborhood will then be 0(exp(— ce^^^'^)) for some c > 0. Thus 
the main contributions to the error come from the remainder terms in the Taylor 
expansions, and this leads to 

Lemma 4.1. [6J Let S^n E C°°{R'^) be a real-valued function, and A^^ E C^{R^). 
Then 



(4.5) 
(4.6) 



|m(0,-)-w^(0,-)|| 



TO ^ £2 



e|hi(0,-)-«^(0,-)|lHi 



< 



Remark 4.1. We note that a cutoff function is necessary and important when one is 
building beams of higher accuracy 

We now show the initial error of time derivative of the GB superposition is also 
under control. We compute the time derivative of fl4.4p to obtain 

dtu'{t,y) = Z{n,e) / dtUGB{t,y; Xo)dxo, 

Jl{0) 

where UGsit, y, xq) = u^^it, y; xq) + W^^it, y; xq) with 



e 



,i*±(i,y;xo)/e 



dtUGBit,y;xo) 

Note that the GB construction ensures that 

dt<^^{t, y; xo) = Tc{y)\V<l>^{t, y; Xo)\ + 0{\y - x(t, Xo)f). 
Recall (ESD we have dtA{t; xo) ~ 0(1). Hence from fl2:T0D we have 

dtUGBiO, y- xo) = [0(1) + e-\B^-'\xo) + 0{\y - xo\))] e^*^(°'^^^«)/^ 



Note that 

Z(n, e) 



HO) L 



0(1) + 



\y - xo\ 



< o(i 



'1/2^ 



which together with Lemma 14.11 again gives 
(4.7) 

||m^(0, ■) - u{0, OIU < elk^(0, ■) - u{0, OIIhi + eWdtU^O, ■) - 9*^(0, OlUg < e'/\ 

Remark 4.2. The above analysis shows that one could choose B^^ to simplify the 
superposition, as was pointed out in Section 2. For example. 



fi) for B 



(4.8) 



(-1) 



-2c(x)|V5in|, then A+{0;Xo) = A^J{xo), A-{0;Xo) = 



u'{t,y) = Zin,e) 



u 



ii) for b[~^^ = 0, then A^{0; Xq 



n+(t) 
ML^^o) 



PGB 



6{w+)dX 



(4.9) 



u'{t,y) = Z{n,e) 



n+{t) 



UpQ^6{w^)dX + / UpQ^6{w )dX 



n-{t) 
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5. Propagation of the approximation error 

We now turn to quantify the evolution error P[m'^]. Recall the Schur's lemma: If 
[Tf]iy) = J K{x,y)f{x)dx and 



then 



sup^ / \K{x,y)\dy = Ci, supy J \K{x,y)\dx = C2, 

IIT/IU. < V^II/IU^. 



Proof. We have by Schwartz 

\[Tf]{y)\'< if \K{x,y)\f{x)dx\ < I \K{x,y)\dx I \K{x,y)\\f{x)\'dx 



<C2j \K{x,y)\\f{x)\'dx. 

So integrating both sides in y and taking the square root gives the result. □ 
We now apply Schur's lemma to a typical term in Jj^^^ P[u'^]dxQ: 

[TA]iy)= [ A{t-x,)F{t,y-xo)e'^^''y-^^'^'^dxo, 

Jl{0) 

where the imaginary part of y; Xq) is bounded below by cl and for convenience 
we will assume that |F| < \y - x{t,xo)\''. Then one can apply Schur's lemma with 



Cl = sup^^ / \y - x(t,Xo)|''e-('=/^)IJ'-^(*'^°)l dy = e^+t / l^l'^e"^!"! dz, and 
(5.1) C2{t,e) = snpy [ |y - xo)| V(^/^)l^-"(*'"°)''rfxo. 



In general one does not know what C2(t,e) will be. As long as A has compact 
support C2 will be at least bounded by ce^^^. Thus the error in norm will be 
bounded by ce'^/^^"''^. We now show that for the wave equation, a better rate can 
be obtained. 

Lemma 5.1. We have 

Proof. From fl2.7p and taking Pq = VxSin{xo) it follows 

^(t,x(t,xo)) = Si„(xo), Vt>0. 
Differentiation of this equation in xq gives 

dx 

-p = po, pit, Xq) ■■= VxS{t, x{t, Xq)). 



dxo 

For non-constant initial phase, at least one element in the deformation matrix ^ 
is non-zero. Assume 7^ near Xq, then writing xq = {xqi^Xq) there exists a 
function h such that Xqi = h{t, z, Xq) and 

z = Xi{t,h{t,z,xo),Xo) 
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in the neighborhood of Xq. Also the map (xqi = h{t, z, xo),xq) {z, xq) is invertible, 
with the Jacobian determined by 



J = det 



/ d{xoi,Xo)\ _ 


dh 




dxi 


V d{z,xo) ) 


dz 




dxQi 



With this map we rewrite the underlying quantity as 

C2 = / (|y-x(t, 2;,xo)P+|i/i-2;n^'/^exp ( - x(t, xo)|^ + \yi - z\^)) Jdxodz. 

Using a stretched coordinate in z so that z — yi = i/e^, with a := y — x(t, z, xq), we 
obtain 

C2 = V~e [ (|a|2 + e|er)'=/V^I«l%xp(--|a|2) Jrfxorfe 

Rewriting e~'^'^'^ = e~'^'^'^/^ ■ e^'^'^'^^^, and using the fact that e"*^'^'^/^ < 1 and 
|^|2g-c|c|2/2 < we obtain 

C2<yfe I (|a|2 + Ce)'^/2g-c|cp/2g-c|apA j^^^^^ 



Here we have used the fact that (|ap + Ce)''/^e"''l"l^/'' < e''/^ for any a e M"~^. As 
long as the initial domain for xq is finitely compact, the above integral is uniformly 
bounded. Note that the local feature of the used map is not restricted, since one 
could use a partition of unity to decompose C2 into a finite sum of terms with the 
same rate of error. The desired estimate thus follows. □ 

This lemma enables us to conclude the following key estimate 
(5.2) ||T[A]|U2 <e^-/2+(i+")/4, 



which will be used to prove the following theorem. 



Theorem 5.2. Let P = — c^{x)A be the linear wave operator and be defined 
m (g;^ with Im{M^) = pi and Z{n,e) = (/5/(27re))"/2. // both An and An have 
compact supports, then u" is an asymptotic solution and satisfies 



(5.3) 



\\P[u']{t,-)\\Ll<e~'^. 
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Proof. Using the volume-preserving map of X = X(t,Xo) and w(t, Xq)) 
w{0,Xo), we obtain 

u'it, y) = Zin, e) / upcBit, y, X(t, Xo))5(w;(t, X{t, Xo)))dXo 
Jn{o) 



Z{n,e) / UGB{t,y;Xo)6{w{0,Xo))dXo 
Jn{o) 

Z{n, e) / UGB{t, y, Xo)6{po - Vr,Sin{xo))dXo 
Jn(o) 



= Z{n,e) UGBit,y;xo)dxo. 
7/(0) 

According to the GB construction, UGB{t,y] xq) is an asymptotic solution for each 
Xo, so will be their superpositions u'^{t,y). It remains to verify (15.31) . First we see 
that 

P[u'{t,y)] = Z{n,e) / P[uGB{t,y] Xo)]dxo, 

7/(0) 

where 

(5.4) P[A{t; xo)e'*(*'^'"'')/^] = {e-\.2{t, y) + e-^c_i + cq) e^*(*'^'-o)/^ 
where for G(t,y) = |9j<l>p - c^|Vy$p, we have 

c^2{t,y) = -G{t,y)A, 



c-iit,y) = 2i 

co{t,y) = dfA{t; xq). 
Using Taylor expansion around x = x(t, Xq) we have 

G{t, y) = G{t, x) + dMt, x)-{y-x) + ^{y- xf dlG{y -x) + 0{\y- a;^. 

Then the Gaussian beam construction sketched in Section 2 ensures that 

\c.2{t,y)\<G\A\\y-x\\ 
Also using the construction for A , we are able to show 

\c-i{t,y)\<G\A\\y-xl \co{t,y)\ < G\A\. 
The construction with positive Im{M) guarantees that 

^it,y;xo) > c\y - x\^. 

Consequently, 



Z-'\\P[u^{t,.)]U.< 



7/(0) 



LI 







i=-2 



\A\\cAe-'^y-''^'''"'^\'^'dxr 



1(0) 



LI 
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continuing the estimate by using the key estimate (15.21) with k = 3,1,0 for F = 
c_2, c_i, Co, respectively 

< ^ ^-1 . ^1/2 ^ ^ ^(l+n)/4 

< -l/2+(l+n)/4 

which when using Z ~ e""-/^ proves the result. □ 

This combined with the obtained initial error and total error estimate in Lemma 
3.1 gives 

Theorem 5.3. Given T > 0, and let u he the solution of the wave equation subject 
to the initial data (u,Ut){0, x) = {A^, Bia)e^^'"^^^^^'^ . Let u'^ he the first order approx- 
imation defined in ^4-^ ''^'i'th initial data satisfying 5*^(0; x) = Sij^{x), M=''(0;x) = 

dlSUx)+ipI, andA^{0-x) = \ ± with \supp{Al^)\ + \supp{Bl^)\ < 

oo. Then there exists eo > 0, a normalization parameter Z(n,e) = ^ , and a 

constant C such that for all e G (0, eo) 



1 I l-n 
4 



(5.5) \\{u^-u){t,-)\\E<Ce-^- 
/orte [0,T]. 

6. An EXAMPLE 

Consider the initial value problem in M'^ for d^u — Au = with initial data 
u(0,x) = e^l^l/^4^> and Ut{0,x) = 0, 

where f{s) G C^(0, oo). Setting g{s) = /(s) exp(zs/e) for s > 0, we extend g{s) to 
be odd on M, i.e. 

9{-\x\) = -f{\x\)exp{i\x\/e). 
This problem has the exact solution 

u{t, x) = 7^{g{t + \x\) - g{t - \x\)). 
\x\ 

At X = this solution has a caustic of the maximum possible strength, since all rays 
starting inward from the sphere \x\ = a arrive at a; = when t = a. This is reflected 
in the behavior of the exact solution 

u{0,t) = g'it) = iifit)/e + f'{t))e^'/% 

which grows like as e goes to zero. 

To build a Gaussian beam approximation for this we need 

ucB{t,x) = \ (^^^''' j^A-^{t,y)e'''^^^^^^^ 
where A^{Q,y) = fi\y\)/\y\ and 

$^(0, x; y) = \y\ + (x - y) ■ p{y) + (x - y) • (^^(J - Piy) + {x - y)/2, 
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where p{y) = y/\y\ and P{y) is the orthogonal projection on the span of p{y). We 
also want $^(0, x; y) = — $^(0, x; y), so that dtu{0, x) = and y) exp(ifc$^(t, x; y)) 
must be a lowest order Gaussian beams concentrated on the null bi-characteristics 
for r ± 1^1 . With these definitions we have 

X] y) = \y\ + {x- x^{t, y)) ■ p{y) + 

\{x - xHt, y)) ■ miy) + ^^T^TT^^^ " ^^^^^^^"^ " 

where a;^(t, x;y)=y± tp{y). For the amplitudes we have 

A^{t, X ± tp{y)) = (1 ± t(l + iP)y^A{0; x). 

Evaluating UGB{t,x) analytically looks difficult, but for ugb(^,0) one has for t > 

UGBit,0)=uit,0) + o{l/e). 

This behavior is predicted by the basic result that, like Fourier integral operators, 
Gaussian beam superpositions give accurate leading order terms in asymptotic ex- 
pansions. 

In principle, one can evaluate the Gaussian beam superposition and compare it 
with the exact solution. Doing this numerically could lead to interesting results on 
the accuracy of these superpositions. 



7. Higher order Approximations 



The accuracy of the phase space based Gaussian beam superposition depends on 
accuracy of the individual Gaussian beam Ansatz. Gaussian beams can be con- 
structed to satisfy the wave equation modulo errors of order e^, for arbitrary N, by 
computing higher order terms in the spatial Taylor series for the phase and ampli- 
tude about the central ray. If we refer the construction in previous sections as the 
first order GB solution, then a k^^ order GB solution will include the Taylor series 
up to {k + iy^ order for the phase, and {k — l — 2iy^ order for the amplitude Ai for 
/ = 0, ■ ■ ■ , [^y^] . The equations for these phase and amplitude Taylor coefficients 
are derived recursively, starting with the phase and then progressing through the 
amplitudes. At each stage (phase function, leading amplitude, next amplitude ...) 
one has to derive the Taylor series up to sufficiently high order before passing to the 
next function in the expansion. 

Let X = X'^it] Xq), with x = x^{t; Xq), denote the bicharacteristic at time t > 0, 
which originates from Xq. Following [6[ [3] we define the k^^ order Gaussian beams 
as follows 



(7.1) 



p{y-x) 



1=0 



exp -T,-_,iM(y) 



where [/](?/) is the k^^ order Taylor polynomial of / about x evaluated at y, and 
p is a cut-off function such that on its support the Taylor expansion of $^ still has 
a positive imaginary part. 
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By invoking the volume preserving map X = X^{t,Xo) and its mverse map 
denoted by Xq = X^{t, X), we obtain a phase space based /c*'' order Gaussian beam 
Ansatz 

Proceeding as previously, we form the superpositions. 



(7.2) ul{t,y) = Z{n,e) 



Jn+(t) Jn-{t) 



where Q{t) = X{t, f2(0)), and X) is the solution of the Liouville equation with 
H = ±c{x)\p\ subject to i(7^(0, X) = p — VxSin{x). 

In (17.11) the initial data for the amplitudes must be chosen consistently with 
the initial data (II. 2p . This leads to the recursion relations 

(7.3) Ai{0,x) + Ar{0,x) = A^ 

dtAt,{0,x) + dtA-^,iO,x)-t{Al{0,x)-A;{0,x))c{x)\VSUx)\ = Bt'\x). 

Note that, since this recursion involves the initial time derivatives of the amplitudes, 
it becomes quite complicated as / increases. 

This gives a k*^^ order asymptotic solution of the wave equation. More precisely, 
we have the following theorem. 

Theorem 7.1. Let P be the linear wave operator of the form P = df — c^A, and u'^ 
IS defined m (T^ with Im{M^) = (31 and Z{n,e) = (/?/(27re))"/^ /3 > 0, then 
is an asymptotic solution and satisfies 

(7.4) lim](i,-)IUg<6^^+V. 

Proof. For notational convenience we estimate only one of two Gaussian beams with 
± index omitted: 



ul{t,y) = Z{n,e) / UkGB{t,y; Xo)dxo. 
Jim 



(0) 

According to the GB construction, UkGB(t,y; xq) are asymptotic solutions for each 
xo, so will be their superpositions ul{t,y). It remains to verify (17.41) . First we see 
that 

Piulit, y)] = Z{n, e) / P[ukGB{t, y; xo)]dxo. 
Ji(o) 



1(0) 

LVJ 



Using dO]) in Section 2 with A replaced by p{y - x) J^Lo ^^T^-i-2i[Miy) 
$ by T^_^_^[^]{y), we have 



and 



c.2it,y) = -Gpiy-x)T^^,[Ao]iy): 

c.,{t,y) = 2iL [pT,-„JAo](t/)] + Gr,%[A](i/), 

ci{t,y) =2iL [pT^3_2,[A,+i](y)]] + Gpr,-_5_2/[A+2](y) + P[pT^t-i-2i[M{y)l ^ = 0, 1, 
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where G = mT^^Miv)? " c^(V,T,-+J$](2/))^]. Using T-_,,[$](y) =^$(y) + 
Rk+i[^]{y)j here -R^+i denotes the remainder of the Taylor expansion, and G{t,y) = 
0{\y — xl'^'^^) we can see that 

\c.2{t,y)\<C\y-x\'+'. 
Also using the construction for Ai and their derivatives, we are able to show 

Ht,y)\<C\y-x\'-'-'^, 

where we have used the fact that differentiation of p vanishes in a neighborhood 
of X. The use of the cut-off function ensures that we can always choose a small 
neighborhood of x(t,xo) so that 



ImiT^^Miy)) > c\y 



x\ 



|2 



Consequently, 

Z-'\\P[u^{t,-)]\\L2 < 



Jl(0) 



dxr 



J=-2 



(0) 

■2 



Li 



LI 



'I(O) 

continuing the estimate by using the key estimate (15.21) 

< ^g-2^fe/2+l _^ ^-1 . ^ j g(l+n)/4 

< fc/2-l+(l+n)/4 

which when using Z ~ e~"/2 proves the result. □ 

In order to obtain an estimate of — u){t,-)\\E for any < t < T, all that 
remains to verify is that the superposition (14. 2 p accurately approximates the initial 
data. However, using the recursion (17.31) to determine the amplitudes, this is again 
an application of [B]. which shows that the initial error in energy norm is bounded 
by e'^/^ for k > 1. Thus our main result for fc*'^ order phase space GB superposition 
is as follows. 

Theorem 7.2. Given T > 0, and let u he the solution of the wave equation sub- 
ject to the initial data {u,ut){0, x) = (Am, -Bin)e*'^'"*^^-'/^, and be the fc*^ order 
approximation defined in Iji7.2\ ) with initial data chosen as described above with 
\supp{Aij^)\ + \supp{Bin)\ < oo. Then there exists eo > 0, a normalization parameter 
Z(n, e) ~ e^"/^, and a constant G such that for all e E (0, cq) 



l-n 
4 



Uu^-u){t,.)\\E<Ge^- 
forte[0,T]. 
Remarks 

• Due to the property of symmetry in time, all results obtained apply to |t| < 
T. 

• For higher order constructions, the Liouville equation for higher order GB 
components can be given similarly to those for the first order GB method. 
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• For computation of high order derivatives of the phase through level set 
functions we refer to [3] for details. 
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